//Purpose: This file does the analysis with the Susenas data;
//This version: September 21, 2020; Jennifer Seager

//NOTE THAT THE SUSENAS DATA IS NOT PUBLICLY AVAILABLE AND CAN BE PURCHASED
//FROM THE INDONESIAN STATISTICAL AGENCY AT https://www.bps.go.id/

//This do file starts from the data file that we generated from the SUSENAS 
//and shows the analysis we did to generate Figrue IV and Table I.II.

//The following is a list of variables that we generate from the Susenas.

/*
year 		year of the data (2010-2013)
female 		=1 if female
age 		age of respondent
malang 		=1 if Malang/ =1 if Pasuruan/Batu
condom 		=1 if use condom with partner (among married women)
work 		=1 if work in last week
income (inv_inc) earnings last month from main job (Inverse-hyperbolic sin transformation)
ltprimary 	=1 if less than primary education
other 		=1 if experienced symptom in last month
noed 		=1 if no education
nosec 		=1 if no secondary education
edyrs 		years of education
child 		=1 if have child (among married women)
married 	=1 if married
nevermarried =1 if never married
divwid		=1 if divorced or widowed
*/

#delimit;
clear;
clear matrix;
eststo clear;
estimates clear;
set more off;
cap log close;

***set directory in file 0;


#delimit cr
/************/
/*Figure IV*/
/************/
use "$da/Susenas-Constructed_2010-2013.dta", clear


keep if age>=28 & age<=42
//keep if age>19 & age < 51
drop if edyrs > 12
keep if female==1

//recreating sample;
gen count =1
egen sampleyear = sum(count),by(year)

set seed 200512
bysort year edyrs (weight cluster female age malang): gen edrand = runiform()

bysort year edyrs (edrand): gen order = sum(count)

gen insample = .

local 0 = "0.0653"
local 1 = "0.0594"
local 2 = "0.0535"
local 3 = "0.0515"
local 4 = "0.0574"
local 5 = "0.0594"
local 6 = "0.3366"
local 7 = "0.0515"
local 8 = "0.0297"
local 9 = "0.1485"
local 10 = "0.0198"
local 11 = "0.0119"
local 12 = "0.0554"

foreach yr in 2010 2011 2012 2013 {
forvalues i = 1/12 {	
	gen temp = sampleyear*``i'' if year==`yr' & edyrs==`i'
	replace insample = 1 if order<=temp & year==`yr' & edyrs==`i'
	drop temp
}
}

//creating age bins


gen agegroup = 1 if age <= 25
replace agegroup = 2 if age >25 & age <=30
replace agegroup = 3 if age >30 & age <=35
replace agegroup = 4 if age >35 & age <=40
replace agegroup = 5 if age >40 & age <=45
replace agegroup = 6 if age >45 & age <=50

gen edgroup = 1 if edyrs<=6
replace edgroup = 3 if edyrs>6 & edyrs<=9
replace edgroup = 4 if edyrs>9 & edyrs<=12

egen strata = group(agegroup edgroup)

local 1 "0.0491"
local 2 "0.0424"
local 3 "0.0089"
local 4 "0.1038"
local 5 "0.0614"
local 6 "0.0134"
local 7 "0.1719"
local 8 "0.0502"
local 9 "0.0290"
local 10 "0.1350"
local 11 "0.0424"
local 12 "0.0257"
local 13 "0.1317"
local 14 "0.0335"
local 15 "0.0045"
local 16 "0.0882"
local 17 "0.0067"
local 18 "0.0022"

gen insample2 = .
foreach yr in 2010 2011 2012 2013 {
forvalues i = 1/18 {	
	gen temp = sampleyear*``i'' if year==`yr' & strata==`i'
	replace insample2 = 1 if order<=temp & year==`yr' & strata==`i'
	drop temp
}
}

keep if insample2==1
 //keep if insample==1

 save "$dt/temp.dta", replace

 
foreach var in  work inv_inc condom other {
reg `var' m_2010 m_2011 m_2012 malang yr2010 yr2011 yr2012 , robust
parmest, saving("$dt/`var'_trends_fem_pb.dta", replace)
}

#delimit;
#delimit;
foreach i in pb {;
use "$dt/work_trends_fem_`i'.dta", clear;
gen reg = "work_fem_`i'";

append using "$dt/condom_trends_fem_`i'.dta";
replace reg = "condom_fem_`i'" if reg=="";

append using "$dt/inv_inc_trends_fem_`i'.dta";
replace reg = "inv_inc_fem_`i'" if reg=="";

append using "$dt/other_trends_fem_`i'.dta";
replace reg = "other_fem_`i'" if reg=="";

gen year = 2010 if parm=="m_2010";
replace year = 2011 if parm=="m_2011";
replace year = 2012 if parm=="m_2012";
replace year = 2013 if parm=="m_2013";
drop if year==.;

save "$dt/eventstudy_data_`i'.dta", replace;
};

use "$dt/eventstudy_data_pb.dta", clear;

twoway scatter estimate year if reg=="condom_fem_pb", mcolor(black)||
	   rspike min95 max95 year if reg=="condom_fem_pb",lpattern(dash) lcolor(black) 
	   || scatteri 0 2013, mcolor(black) msymbol(circle) msize(medium)
	   sort(year) yline(0, lcolor(gs12)) legend(off) 
	   ytitle("Coefficient on Malang*Year") xtitle("Year") title("Condom Use")
	   scheme(s1mono) graphregion(color(white)) plotregion(lwidth(vvthin))
	   ylabel(-.1(0.02)0.1) ymtick(-.1(0.01).1)
	   xlabel(2010(1)2013) xmtick(2010(1)2013);
graph export "$dout/Fig_IV_condomuse.pdf", as(pdf) replace;

twoway scatter estimate year if reg=="work_fem_pb", mcolor(black)||
	   rspike min95 max95 year if reg=="work_fem_pb",lpattern(dash) lcolor(black) 
	   || scatteri 0 2013, mcolor(black) msymbol(circle) msize(medium)
	   sort(year) yline(0, lcolor(gs12)) legend(off) 
	   ytitle("Coefficient on Malang*Year") xtitle("Year") 
	   title("Worked in last week")
	   scheme(s1mono) graphregion(color(white)) plotregion(lwidth(vvthin)) 
	   ylabel(-1(0.2)1) ymtick(-1(0.1)1)
	   xlabel(2010(1)2013) xmtick(2010(1)2013);
graph export "$dout/Fig_IV_worklastweek.pdf", as(pdf) replace;

twoway scatter estimate year if reg=="other_fem_pb", mcolor(black)||
	   rspike min95 max95 year if reg=="other_fem_pb",lpattern(dash) lcolor(black) 
	   || scatteri 0 2013, mcolor(black) msymbol(circle) msize(medium)
	   sort(year) yline(0, lcolor(gs12)) legend(off) 
	   ytitle("Coefficient on Malang*Year") xtitle("Year") 
	   title("Experienced Health Symptom in Past Month")
	   scheme(s1mono) graphregion(color(white)) plotregion(lwidth(vvthin)) 
	   ylabel(-1(0.2)1) ymtick(-1(0.1)1)
	   xlabel(2010(1)2013) xmtick(2010(1)2013);
graph export "$dout/Fig_IV_othersymptomlastmonth.pdf", as(pdf) replace;

twoway scatter estimate year if reg=="inv_inc_fem_pb", mcolor(black)||
	   rspike min95 max95 year if reg=="inv_inc_fem_pb",lpattern(dash) lcolor(black) 
	   || scatteri 0 2013, mcolor(black) msymbol(circle) msize(medium)
	   sort(year) yline(0, lcolor(gs12)) legend(off) 
	   ytitle("Coefficient on Malang*Year") xtitle("Year") 
	   title("Income earned in last month")
	   scheme(s1mono) graphregion(color(white)) plotregion(lwidth(vvthin)) 
	   ylabel(-10(2)10) ymtick(-10(1)10)
	   xlabel(2010(1)2013) xmtick(2010(1)2013);
graph export "$dout/Fig_IV_inv_inc.pdf", as(pdf) replace;


/************/
/*Table I.II*/
/************/;
#delimit;
use "$da/Susenas-Constructed_2010-2013.dta", clear;

keep if age>=28 & age <=42;
keep if year==2013;
preserve;
use "$da/Crimes-Against-Morality_Main_data-constructed.dta", clear;

keep if dataset==1& wave==1 & married !=.;
keep if age>=28 & age<=42;
gen count =1;
egen agenum = sum(count), by(age);
egen total = sum(count);

gen agepc = agenum/total;
tab agepc ; 

duplicates drop age, force;
keep age agenum total agepc count;

save "$dt/agepctile.dta", replace;

restore;

merge m:1 age using "$dt/agepctile.dta";
keep if female==1;
keep if _merge==3;

egen umurnum = sum(count),by(age);
egen total2  = sum(count);
gen umurpc = umurnum/total2;

gen pcdiff = umurpc - agepc;
gen pcdiff1 = abs(pcdiff);

set seed 180410;
bysort age (cluster female weight income ltprimary): gen rand = runiform();

sort age rand;

bysort age (rand): gen id = sum(count);

gen cut = round(agepc*total2*.55);
	//multiplying by .55 becuase this makes it so that we are pulling a subset
	//of women from every age group so that we can acheive the desired age distribution;

gen keep = 1 if id<=cut;
replace keep=0 if keep==.;

***For Table I.II (column 2);
mean married divwid nevermarried edyrs age [pw=weight] if female==1 &keep==1;

mean child [pw=weight] if female==1 &keep==1;

use "$da/Crimes-Against-Morality_Main_data-constructed.dta", clear;


***For Table I.II (column 1);
mean married divwid unmarried edyrs age child if wave==1;


